This notebook pulls in data on Nipah RBP entry into CHO-EFNB2 and CHO-EFNB3 cells, filters data, calculates stats, and makes figures¶
In [1]:
# this cell is tagged as parameters for `papermill` parameterization
e2_distances_file = None
func_scores_E2_file = None
func_scores_E3_file = None
filtered_E2_data = None
filtered_E3_data = None
contact_type_plot = None
E2_E3_entry_corr_plot = None
E2_E3_entry_all_muts_plot = None
combined_E2_E3_correlation_plots = None
entry_region_boxplot_plot = None
nipah_config = None
altair_config = None
entropy_file = None
surface = None
In [2]:
# Parameters
func_scores_E2_file = "results/func_effects/averages/CHO_EFNB2_low_func_effects.csv"
func_scores_E3_file = "results/func_effects/averages/CHO_EFNB3_low_func_effects.csv"
e2_distances_file = "results/distances/2vsm_distances.csv"
contact_type_plot = "results/images/contact_type_plot.html"
filtered_E2_data = "results/filtered_data/E2_entry_filtered.csv"
filtered_E3_data = "results/filtered_data/E3_entry_filtered.csv"
E2_E3_entry_corr_plot = "results/images/E2_E3_entry_corr_plot.html"
E2_E3_entry_all_muts_plot = "results/images/E2_E3_entry_all_muts_plot.html"
combined_E2_E3_correlation_plots = (
"results/images/combined_E2_E3_correlation_plots.html"
)
nipah_config = "nipah_config.yaml"
altair_config = "data/custom_analyses_data/theme.py"
entropy_file = "results/entropy/entropy.csv"
entry_region_boxplot_plot = "results/images/entry_region_boxplot_plot.html"
surface = "data/custom_analyses_data/surface_exposure.csv"
In [3]:
import math
import os
import re
import altair as alt
import numpy as np
import pandas as pd
import scipy.stats
import Bio.SeqIO
import yaml
from Bio import AlignIO
from Bio import PDB
from Bio.Align import PairwiseAligner
In [4]:
# allow more rows for Altair
_ = alt.data_transformers.disable_max_rows()
if os.getcwd() == '/fh/fast/bloom_j/computational_notebooks/blarsen/2023/Nipah_Malaysia_RBP_DMS/':
pass
print("Already in correct directory")
else:
os.chdir("/fh/fast/bloom_j/computational_notebooks/blarsen/2023/Nipah_Malaysia_RBP_DMS/")
print("Setup in correct directory")
Setup in correct directory
In [5]:
if surface is None:
e2_distances_file = "results/distances/2vsm_distances.csv"
func_scores_E2_file = "results/func_effects/averages/CHO_EFNB2_low_func_effects.csv"
func_scores_E3_file = "results/func_effects/averages/CHO_EFNB3_low_func_effects.csv"
filtered_E2_data = "results/filtered_data/E2_entry_filtered.csv"
filtered_E3_data = "results/filtered_data/E3_entry_filtered.csv"
nipah_config = "nipah_config.yaml"
altair_config = "data/custom_analyses_data/theme.py"
entropy_file = "results/entropy/entropy.csv"
surface = "data/custom_analyses_data/surface_exposure.csv"
Read in custom altair theme and Import YAML file with parameters¶
In [6]:
if altair_config:
with open(altair_config, 'r') as file:
exec(file.read())
with open(nipah_config) as f:
config = yaml.safe_load(f)
Filter and merge EFNB2 and EFNB3 entry¶
In [7]:
#Import median entry scores from different selections
func_scores_E2 = pd.read_csv(func_scores_E2_file).dropna().round(3) # this removes wildtype observations from df
func_scores_E3 = pd.read_csv(func_scores_E3_file).dropna().round(3) # this removes wildtype observations from df
In [8]:
def num_selections_and_filter(df, name):
def calculate_filtering_mutants(df, name):
print(f"\nThe dataset is: {name}")
# How many selections were performed
max_sels = df["n_selections"].max()
num_sels_cutoff = (max_sels / 2) + 1
# Find total number of possible mutants
total_mut = (602 - 71) * 19
# Filter data
filter_test = df[
(df["site"] != 603) & (df["mutant"] != "-") & (df["mutant"] != "*")
]
num_variants_pre_filter = filter_test.shape[0]
print(
f"After filtering stop and gaps, there are {num_variants_pre_filter} mutants which is {(num_variants_pre_filter/total_mut) * 100:.1f}%"
)
filter_test_times_seen = filter_test[
filter_test["times_seen"] >= config["func_times_seen_cutoff"]
]
num_variants_times_seen = filter_test_times_seen.shape[0]
print(
f"After filtering for {config['func_times_seen_cutoff']} times seen, there are {num_variants_times_seen}, which is {(num_variants_times_seen/total_mut) * 100:.1f}%"
)
filter_test_effect_std = filter_test_times_seen[
filter_test_times_seen["effect_std"] <= config["func_std_cutoff"]
]
num_variants_std = filter_test_effect_std.shape[0]
print(
f"After filtering for {config['func_std_cutoff']} std cutoff, there are {num_variants_std}, which is {(num_variants_std/total_mut) * 100 :.1f}%"
)
filter_test_n_selections = filter_test_effect_std[
filter_test_effect_std["n_selections"] >= num_sels_cutoff
]
num_variants_n_selections = filter_test_n_selections.shape[0]
print(
f"After filtering for mutants in in all selections, there are {num_variants_n_selections}, which is {(num_variants_n_selections/total_mut) * 100 :.1f}%"
)
def apply_filters(df, name):
# Now do the main filtering
max_sels = df["n_selections"].max()
num_sels_cutoff = (max_sels / 2) + 1
print(
f"The number of selections a mutant must be observed in is: {num_sels_cutoff}"
)
# The main filtering. Filters site 603 (is a stop codon/end of gene and we don't want those mutants). Also filter out stop mutants and apply filtering from config file
filtered_df = df[
(df["site"] != 603)
& (df["mutant"] != "-")
& (df["mutant"] != "*")
& (df["times_seen"] >= config["func_times_seen_cutoff"])
& (df["effect_std"] <= config["func_std_cutoff"])
& (df["n_selections"] >= num_sels_cutoff)
]
return filtered_df
# Filtering stats
calculate_filtering_mutants(df, name) # call definition above
# Filter
filtered_df = apply_filters(df, name)
# Now write filtered_data to results/
if name == "E2":
filtered_df.to_csv(filtered_E2_data, index=False)
if name == "E3":
filtered_df.to_csv(filtered_E3_data, index=False)
# return filtered dataframe
return filtered_df
# Call the filtering functions
func_scores_E2 = num_selections_and_filter(func_scores_E2, "E2")
func_scores_E3 = num_selections_and_filter(func_scores_E3, "E3")
# make a merged dataframe of ephrin-b2 and ephrin-b3 entry data
def merge_data(df1,df2):
merged_df = pd.merge(
df1,
df2,
on=["site", "mutant", "wildtype"],
how="outer",
suffixes=["_E2", "_E3"],
)
df1["cell_type"] = "CHO-EFNB2"
df2["cell_type"] = "CHO-EFNB3"
concat_df = pd.concat([df1,df2])
return merged_df,concat_df
merged_df,concat_df = merge_data(func_scores_E2,func_scores_E3)
# Show some stats of filtered merged data
stats = merged_df.describe().round(2)
display(stats)
The dataset is: E2 After filtering stop and gaps, there are 10040 mutants which is 99.5% After filtering for 2 times seen, there are 9789, which is 97.0% After filtering for 1 std cutoff, there are 9729, which is 96.4% After filtering for mutants in in all selections, there are 9725, which is 96.4% The number of selections a mutant must be observed in is: 5.0 The dataset is: E3 After filtering stop and gaps, there are 10072 mutants which is 99.8% After filtering for 2 times seen, there are 9852, which is 97.7% After filtering for 1 std cutoff, there are 9714, which is 96.3% After filtering for mutants in in all selections, there are 9586, which is 95.0% The number of selections a mutant must be observed in is: 4.5
| site | effect_E2 | effect_std_E2 | times_seen_E2 | n_selections_E2 | effect_E3 | effect_std_E3 | times_seen_E3 | n_selections_E3 | |
|---|---|---|---|---|---|---|---|---|---|
| count | 9885.00 | 9725.00 | 9725.00 | 9725.00 | 9725.00 | 9586.00 | 9586.00 | 9586.00 | 9586.00 |
| mean | 337.22 | -0.92 | 0.32 | 7.20 | 7.99 | -0.95 | 0.34 | 6.14 | 6.98 |
| std | 153.23 | 1.40 | 0.24 | 4.08 | 0.12 | 1.46 | 0.25 | 3.27 | 0.14 |
| min | 71.00 | -3.52 | 0.00 | 2.00 | 5.00 | -3.61 | 0.00 | 2.00 | 5.00 |
| 25% | 205.00 | -2.21 | 0.15 | 4.62 | 8.00 | -2.50 | 0.17 | 4.14 | 7.00 |
| 50% | 338.00 | -0.27 | 0.30 | 6.38 | 8.00 | -0.22 | 0.33 | 5.43 | 7.00 |
| 75% | 470.00 | 0.22 | 0.47 | 8.50 | 8.00 | 0.22 | 0.49 | 7.14 | 7.00 |
| max | 602.00 | 0.64 | 1.00 | 64.38 | 8.00 | 0.66 | 1.00 | 49.14 | 7.00 |
Stats¶
In [9]:
def calculate_stats(df, name):
print(f"For {name}:")
total_mut = (602 - 71) * 19
print(f'There are {total_mut} amino acid mutations possible')
muts_present = df["effect"].shape[0]
fraction_muts = muts_present / total_mut
print(
f"fraction muts present in {name} is {fraction_muts:.2f} {muts_present}/{total_mut}"
)
deleterious_muts = df[df["effect"] <= -0.25].shape[0]
neutral_muts = df[(df["effect"] > -0.25) & (df["effect"] < 0.25)].shape[0]
positive_muts = df[df["effect"] > 0.25].shape[0]
frac_bad_muts = deleterious_muts / muts_present
frac_neutral_muts = neutral_muts / muts_present
frac_pos_muts = positive_muts / muts_present
print(
f"The number of deleterious mutants for {name} is {frac_bad_muts:.2f} {deleterious_muts}/{muts_present}"
)
print(
f"The number of neutral mutants for {name} is {frac_neutral_muts:.2f} {neutral_muts}/{muts_present}"
)
print(
f"The number of positive mutants for {name} is {frac_pos_muts:.2f} {positive_muts}/{muts_present}\n"
)
calculate_stats(func_scores_E2, "E2")
calculate_stats(func_scores_E3, "E3")
For E2: There are 10089 amino acid mutations possible fraction muts present in E2 is 0.96 9725/10089 The number of deleterious mutants for E2 is 0.50 4907/9725 The number of neutral mutants for E2 is 0.26 2571/9725 The number of positive mutants for E2 is 0.23 2239/9725 For E3: There are 10089 amino acid mutations possible fraction muts present in E3 is 0.95 9586/10089 The number of deleterious mutants for E3 is 0.49 4722/9586 The number of neutral mutants for E3 is 0.28 2680/9586 The number of positive mutants for E3 is 0.23 2178/9586
How many sites and which sites only have negative entry scores for mutations?¶
In [10]:
def overall_stats_all_neg(df,effect):
filtered_df = df.groupby('site').filter(lambda group: (group[effect] < 0).all())
unique = filtered_df['site'].unique()
print(list(unique))
total_sites = df['site'].unique().shape[0]
subset = filtered_df['site'].unique().shape[0]
fraction = subset/total_sites
percent = fraction * 100
print(f'The total number of sites are: {total_sites}')
print(f' The number of sites where all mutants are negative for {effect}: {subset}')
print(f' The percent of sites where all mutants are negative for {effect}: {percent:.0f}')
return unique
intolerant_sites_E2 = list(overall_stats_all_neg(func_scores_E2,'effect'))
intolerant_sites_E3 = list(overall_stats_all_neg(func_scores_E3,'effect'))
[80, 106, 107, 108, 111, 112, 113, 120, 121, 125, 126, 127, 130, 138, 146, 151, 157, 158, 159, 162, 163, 165, 167, 172, 189, 203, 205, 206, 207, 208, 216, 229, 240, 246, 251, 253, 254, 257, 258, 259, 260, 262, 263, 264, 266, 267, 303, 323, 331, 347, 355, 382, 387, 395, 412, 460, 467, 487, 489, 493, 499, 500, 503, 506, 537, 563, 565, 574, 594, 598] The total number of sites are: 532 The number of sites where all mutants are negative for effect: 70 The percent of sites where all mutants are negative for effect: 13 [95, 100, 106, 107, 108, 111, 112, 113, 121, 125, 126, 138, 146, 158, 162, 163, 201, 203, 206, 207, 216, 229, 240, 243, 248, 251, 253, 257, 258, 259, 260, 263, 266, 303, 347, 352, 355, 368, 382, 387, 389, 395, 412, 439, 458, 460, 467, 486, 487, 489, 493, 494, 497, 499, 500, 501, 503, 504, 505, 506, 510, 526, 531, 532, 533, 537, 556, 557, 563, 565, 573, 574, 579, 581, 584, 585, 588, 594] The total number of sites are: 532 The number of sites where all mutants are negative for effect: 78 The percent of sites where all mutants are negative for effect: 15
In [11]:
def calculate_top_func_scores(df,effect):
percentile_95_effect_E2 = df[effect].quantile(0.999)
cutoff_E2_df_sites = df[df[effect] > percentile_95_effect_E2]
E2_site_cutoff = cutoff_E2_df_sites['site'].unique()
print(f'The sites with the highest functional scores for {effect} are: {list(E2_site_cutoff)}')
calculate_top_func_scores(func_scores_E2,'effect')
calculate_top_func_scores(func_scores_E3,'effect')
The sites with the highest functional scores for effect are: [280, 305, 407, 463, 468, 501, 552, 556, 584, 599] The sites with the highest functional scores for effect are: [183, 315, 337, 358, 378, 393, 418, 419, 554, 597]
Make bubble plots of receptor contact site type (median values per site)¶
In [12]:
def make_bubbleplot_entry_region(df): # Create a bubble plot using Altair for contact site mutants
barrel_ranges = {
"Hydrophobic": config["hydrophobic"],
"Salt Bridges": config["salt_bridges"],
"Hydrogen Bonds": config["h_bond_total"],
"Contact": config["contact_sites"],
"Overall": list(range(71, 602)),
}
custom_order = [
"Hydrophobic",
"Salt Bridges",
"Hydrogen Bonds",
"Contact",
"Overall",
]
empty_chart = []
variant_selector = alt.selection_point(
on="mouseover", empty=False, fields=["site"], value=1
)
for selection in ["CHO-EFNB2", "CHO-EFNB3"]:
agg_means = []
tmp_df = df[df["cell_type"] == selection]
mean_df = tmp_df.groupby("site")[["effect"]].median().reset_index()
# For each barrel, filter the site_means dataframe to the sites belonging to that barrel and then store the means
for barrel, sites in barrel_ranges.items():
subset = mean_df[mean_df["site"].isin(sites)]
for _, row in subset.iterrows():
agg_means.append(
{"barrel": barrel, "effect": row["effect"], "site": row["site"]}
)
agg_means_df = pd.DataFrame(agg_means)
chart = (
alt.Chart(agg_means_df, title=f"{selection}")
.mark_point()
.encode(
x=alt.X(
"barrel:O",
sort=custom_order,
title='Contact Type',
axis=alt.Axis(labelAngle=-90),
),
y=alt.Y(
"effect",
title="Median Cell Entry",
axis=alt.Axis(grid=True, tickCount=4),
),
xOffset="random:Q",
tooltip=["barrel", "effect", "site"],
size=alt.condition(
variant_selector, alt.value(100), alt.value(25)
),
color=alt.condition(
variant_selector, alt.value("orange"), alt.value("black")
),
opacity=alt.condition(variant_selector, alt.value(1), alt.value(0.3)),
).properties(width=config['width'],height=config['height'])
.transform_calculate(random="sqrt(-1*log(random()))*cos(2*PI*random())")
)
empty_chart.append(chart)
combined_effect_chart = (
alt.hconcat(*empty_chart)
.resolve_scale(y="shared", x="shared", color="independent")
.add_params(variant_selector)
)
return combined_effect_chart
tmp_img = make_bubbleplot_entry_region(concat_df)
tmp_img.display()
if entry_region_boxplot_plot is not None:
tmp_img.save(contact_type_plot)
Make bubble plots depending on amino acid property¶
In [13]:
def make_bubbleplot_wildtype_prop(df):
barrel_ranges = {
"hydrophobic": list(["A", "V", "L", "I", "M"]),
"aromatic": list(["Y", "W", "F"]),
"positive": list(["K", "R", "H"]),
"negative": list(["E", "D"]),
"hydrophilic": list(["S", "T", "N", "Q"]),
"special": list(["C", "P", "G"]),
}
empty_charts = []
variant_selector = alt.selection_point(
on="mouseover", empty=False, nearest=True, fields=["site"], value=1
)
for selection in ["CHO-EFNB2", "CHO-EFNB3"]:
if selection == "CHO-EFNB2":
effect_name = "EFNB2"
else:
effect_name = "EFNB3"
tmp_df = df[df["cell_type"] == selection]
unique_wildtype_df = tmp_df[["site", "wildtype"]].drop_duplicates()
mean_df = tmp_df.groupby("site")[["effect"]].median().reset_index()
mean_df = pd.merge(mean_df, unique_wildtype_df, on="site", how="left")
agg_means = []
# For each barrel, filter the site_means dataframe to the sites belonging to that barrel and then store the means
for barrel, sites in barrel_ranges.items():
subset = mean_df[mean_df["wildtype"].isin(sites)]
for _, row in subset.iterrows():
agg_means.append(
{"wildtype_class": barrel, "effect": row["effect"], "site": row["site"], "wildtype": row["wildtype"]}
)
agg_means_df = pd.DataFrame(agg_means)
chart = (
alt.Chart(agg_means_df, title=f"{selection}")
.mark_point()
.encode(
x=alt.X(
"wildtype_class:O",
title="Wildtype amino acid property",
axis=alt.Axis(labelAngle=-90),
),
y=alt.Y(
"effect",
title=f"Median Cell Entry",
axis=alt.Axis(grid=True, tickCount=4),
),
xOffset="random:Q",
tooltip=["wildtype_class", "effect", "site","wildtype"],
size=alt.condition(variant_selector, alt.value(100), alt.value(25)),
color=alt.condition(
variant_selector, alt.value("orange"), alt.value("black")
),
opacity=alt.condition(variant_selector, alt.value(1), alt.value(0.2)),
).properties(width=config['width'],height=config['height'])
.transform_calculate(random="sqrt(-1*log(random()))*cos(2*PI*random())")
)
empty_charts.append(chart)
combined_effect_chart = (
alt.hconcat(*empty_charts)
.resolve_scale(y="shared", x="shared", color="independent")
.add_params(variant_selector)
)
return combined_effect_chart
wildtype_aa_bubble_img = make_bubbleplot_wildtype_prop(concat_df)
wildtype_aa_bubble_img.display()
Plot correlations between E2 and E3 entry¶
In [14]:
# Import distance data
e2_distances = pd.read_csv(e2_distances_file)
distance_df = pd.merge(
merged_df, e2_distances[["site", "distance"]], on="site", how="left"
)
def determine_distance(df):
# Define the conditions
conditions = [
df["distance"] < 4,
(df["distance"] >= 4) & (df["distance"] <= 8),
df["distance"] > 8,
]
# Define the associated values for the conditions
choices = ["contact", "close", "distant"]
# Apply the conditions and choices to the 'E2_contact' column
df["contact"] = np.select(conditions, choices, default="distant")
return df
distance_df = determine_distance(distance_df)
def median_correlation_plot(df, metric):
aggregation = getattr(df.groupby("site")[["effect_E2", "effect_E3"]], metric)
means = aggregation().reset_index()
slope, intercept, r_value, p_value, std_err = scipy.stats.linregress(
means["effect_E2"], means["effect_E3"]
)
display(r_value.round(2))
means = means.rename(
columns={"effect_E2": f"{metric}_E2", "effect_E3": f"{metric}_E3"}
)
contact_sites = df[["site", "contact", "wildtype"]].drop_duplicates()
df_mean = pd.merge(means, contact_sites, on="site", how="left")
chart = (
alt.Chart(df_mean)
.mark_point()
.encode(
x=alt.X(f"{metric}_E2", title="Summed Cell Entry for EFNB2"),
y=alt.Y(f"{metric}_E3", title="Summed Cell Entry for EFNB3"),
tooltip=["site", "wildtype"],
color=alt.Color(
"contact",
scale=alt.Scale(
domain=["contact", "close", "distant"],
range=["#1f4e79", "#ff7f0e", "gray"],
),
legend=alt.Legend(title="RBP Distance to Receptor"),
),
)
)
min_ = int(df_mean[f"{metric}_E2"].min())
max_ = int(df_mean[f"{metric}_E3"].max())
text = (
alt.Chart({"values": [{"x": min_, "y": max_, "text": f"r = {r_value:.2f}"}]})
.mark_text(
align="left",
baseline="top",
dx=0, # Adjust this for position
dy=0, # Adjust this for position
)
.encode(x=alt.X("x:Q"), y=alt.Y("y:Q"), text="text:N")
)
plot = chart + text
return plot
E2_E3_plot = median_correlation_plot(distance_df, "sum")
E2_E3_plot.display()
if entry_region_boxplot_plot is not None:
E2_E3_plot.save(E2_E3_entry_corr_plot)
def correlation_plot(df):
variant_selector = alt.selection_point(
on="mouseover", empty=False, nearest=True, fields=["contact"], value=1
)
chart = (
alt.Chart(df)
.mark_point(opacity=0.2,size=15)
.encode(
x=alt.X("effect_E2", title="Cell Entry for EFNB2"),
y=alt.Y("effect_E3", title="Cell Entry for EFNB3"),
tooltip=["site", "wildtype", "mutant"],
color=alt.Color(
"contact",
scale=alt.Scale(
domain=["contact", "close", "distant"],
range=["#1f4e79", "#ff7f0e", "gray"],
),
legend=alt.Legend(title="RBP Distance to Receptor"),
),
opacity=alt.condition(variant_selector,alt.value(1),alt.value(0.2)),
).add_params(variant_selector)
)
min_ = int(df["effect_E2"].min())
max_ = int(df["effect_E3"].max())
return chart
tmp_img_corr = correlation_plot(distance_df)
tmp_img_corr.display()
if entry_region_boxplot_plot is not None:
tmp_img_corr.save(E2_E3_entry_all_muts_plot)
if entry_region_boxplot_plot is not None:
(E2_E3_plot | tmp_img_corr).save(combined_E2_E3_correlation_plots)
0.81
Make boxplot showing median entry by RBP region¶
In [15]:
def make_boxplot_entry_region(df):
barrel_ranges = {
"Stalk": list(range(96, 147)),
"Neck": list(range(148, 165)),
"Linker": list(range(166, 177)),
"Head": list(range(178, 602)),
'Receptor Contact': config['contact_sites'],
"Total": list(range(71, 602)),
}
custom_order = ["Stalk", "Neck", "Linker", "Head", "Receptor Contact", "Total"]
empty_charts = []
for selection in ["CHO-EFNB2", "CHO-EFNB3"]:
if selection == "CHO-EFNB2":
effect_name = "EFNB2"
else:
effect_name = "EFNB3"
tmp_df = df[df["cell_type"] == selection]
agg_means = []
# For each barrel, filter the site_means dataframe to the sites belonging to that barrel and then store the means
for barrel, sites in barrel_ranges.items():
subset = tmp_df[tmp_df["site"].isin(sites)]
for _, row in subset.iterrows():
agg_means.append(
{"region": barrel, "effect": row["effect"], "site": row["site"]}
)
agg_means_df = pd.DataFrame(agg_means)
chart = (
alt.Chart(agg_means_df, title=f"{selection}")
.mark_boxplot(color='darkgray')
.encode(
x=alt.X(
"region:O",
sort=custom_order,
title="RBP Region",
axis=alt.Axis(labelAngle=-90),
),
y=alt.Y(
"effect",
title=f"Cell Entry",
axis=alt.Axis(grid=True, tickCount=4),
),
tooltip=["region", "effect", "site"],
).properties(width=config['width'],height=config['height'])
)
empty_charts.append(chart)
combined_effect_chart = alt.hconcat(*empty_charts).resolve_scale(
y="shared", x="shared", color="independent"
)
return combined_effect_chart
entry_region_boxplot = make_boxplot_entry_region(concat_df)
entry_region_boxplot.display()
if entry_region_boxplot_plot is not None:
entry_region_boxplot.save(entry_region_boxplot_plot)
Check for potential neutral/beneficial glycosylation sites¶
In [16]:
def find_potential_glycan_sites(df):
filtered = df[df["wildtype"].isin(["T", "S"])]
matching_sites = []
for index, row in filtered.iterrows():
# Check for existence of site two numbers prior with 'N' mutant and positive effect
prior_rows = df[
(df["site"] == row["site"] - 2) & (df["mutant"] == "N") & (df["effect"] > 0)
]
if not prior_rows.empty:
matching_sites.append(row["site"])
unique_list1 = list(set(matching_sites))
unique_list1 = [x - 2 for x in unique_list1]
filtered = df[df["wildtype"].isin(["N"])]
matching_sites = []
for index, row in filtered.iterrows():
# Check for existence of site two numbers prior with 'N' mutant and positive effect
prior_rows = df[
(df["site"] == row["site"] + 2)
& (df["mutant"].isin(["T", "S"]))
& (df["effect"] > 0)
]
if not prior_rows.empty:
matching_sites.append(row["site"])
unique_list2 = list(set(matching_sites))
unique_list = unique_list1 + unique_list2
unique_list.sort()
print(f"The total number of potential PNLG sites are: {len(unique_list)}")
return unique_list
PNLG = find_potential_glycan_sites(func_scores_E3)
surface_df = pd.read_csv(surface)
surface_df.rename(columns={"exposure_A": "Surface Exposure"}, inplace=True)
PNLG_SURFACE = surface_df[surface_df["site"].isin(PNLG)]
PNLG_SURFACE = list(
PNLG_SURFACE[PNLG_SURFACE["Surface Exposure"] >= 25]["site"].unique()
)
print(f"\nThe surface exposed PNLG sites are: {PNLG_SURFACE}\n")
glycans = config["glycans"]
filtered_PNLG_SURFACE = [num for num in PNLG_SURFACE if num not in glycans]
print(filtered_PNLG_SURFACE)
print(len(filtered_PNLG_SURFACE))
The total number of potential PNLG sites are: 33 The surface exposed PNLG sites are: [159, 180, 191, 192, 275, 288, 306, 311, 326, 378, 383, 403, 417, 423, 473, 478, 518, 543, 554, 570, 600] [180, 191, 192, 275, 288, 311, 326, 383, 403, 423, 473, 478, 518, 543, 554, 570, 600] 17
In [ ]: